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ABSTRACT 

We determine the torque exerted in a steady state by an external potential on 
a three-dimensional gaseous disk at a non-coorbital corotation resonance. Our 
model accounts for the feedback of the torque on the surface density and vorticity 
in the corotation region, and assumes that the disk has a barotropic equation of 
state and a nonzero effective viscosity. The ratio of the torque to the value 
given by the formula of Goldreich & Tremaine depends essentially on a single 
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dimensionless parameter, which quantifies the extent to which the resonance is 
(T) ' saturated. We discuss the implications for the eccentricity evolution of young 
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fcj I 1. Introduction 

CO 

The forcing of a gaseous disk by a gravitational perturber at a resonance can result in a 
strong response and an interchange of energy and angular momentum between the perturber 
and the disk. A corotation resonance arises where the pattern speed of the forcing matches 
the local angular speed of disk material. This occurs in the context of galaxies where a 
stellar bar or spiral arms force motions in the interstellar medium (Binney & Tremaine 
1987). Corotation resonances also arise when satellites orbit within a disk, as occurs with 
young planets in protoplanetary disks, or moons in planetary ring systems (Goldreich & 
Tremaine 1979, 1980; hereafter GT79 and GT80 respectively). 

For planets, two types of corotation resonance need to be distinguished. A coorbital 
corotation resonance occurs where the orbital period of the planet matches the orbital period 
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of disk material. For a planet with a circular orbit, this is the only type of corotation reso- 
nance that can arise. The second type of corotation resonance is non-coorbital and occurs if 
the planet executes an eccentric orbit. The forcing due to a planet with an eccentric orbit can 
be decomposed into a series of rotating components of various strengths and pattern speeds 
(GT80). Bar or spiral galaxies, on the other hand, give rise to non-coorbital resonances. 

The analyses of the coorbital and non-coorbital resonances are different. In the coor- 
bital case, which has recently been considered by Masset (2001, 2002) and Balmforth & 
Korycansky (2001), the forcing is stronger and involves multiple Fourier components. There 
is more likely to be a reduction in density caused by the tendency of the planet to open 
a gap in the disk. In this paper, however, we are concerned with non-coorbital corotation 
resonances, which are critical in determining the eccentricity evolution of planets resulting 
from planet-disk interactions (GT80), and are therefore likely to be important in explain- 
ing the eccentricities of many of the observed extrasolar planets. For a disk in which a 
gap is opened by a planet having a small orbital eccentricity, corotation resonances tend to 
damp the planet's eccentricity, while Lindblad resonances cause it to grow. The effects of 
each type of resonance are strong, but are nearly equal in magnitude, to the extent that 
they nearly cancel. The balance is slightly in favor of eccentricity damping, if the corota- 
tion resonances operate at maximal efficiency (i.e., are unsaturated). The final outcome of 
eccentricity evolution depends on the details. 

The disturbances of the disk caused by a corotation resonance remain localized to the 
corotation region (GT79). They are unable to propagate away from the resonance, as occurs 
in the case of Lindblad resonances. Instead, they act back locally on the disk and may 
change its density and vorticity in such a manner as to reduce (saturate) the corotation 
torque (e.g., Ward 1991). However, the effective turbulent viscosity of the disk could lessen 
the saturation by limiting the back- reaction on the disk (Ward 1992). Full saturation occurs 
when the torque is reduced to zero, but even a small degree of saturation (5%) could change 
the sense of eccentricity evolution in the GT80 model from decay to growth (Goldreich & 
Sari 2002). 

In this paper, we develop a detailed model for the saturation of a non-coorbital corota- 
tion resonance in a viscous accretion disk. We determine how the strength of the corotation 
torque varies with the strength of the tidal forcing and the effective viscosity of the disk. 
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2. Physical picture 

The interaction between the perturber and the disk in the corotation region can be 
considered to involve two distinct aspects. In one aspect, the tidal potential excites an 
evanescent disturbance in the disk, i.e., a disturbance that decays exponentially with dis- 
tance from the resonance over the corotation region. Associated with this disturbance are 
evanescent angular momentum fluxes that are oppositely directed on the two sides of the 
resonance. In the linear, inviscid theory of GT79, there is a jump in the flux at the res- 
onance that corresponds to a torque being injected into the disk exactly at the corotation 
radius. The torque depends on the radial derivative of the ratio of the surface density to the 
vorticity at corotation. In the second aspect, the disk responds to the torque by changing 
the distributions of surface density and angular velocity over the corotation region (Lubow 
1990). These changes are such as to cause a reduction in the torque, and, ultimately, it is 
expected that the torque will be eliminated. This picture is time-dependent because the 
torque injection continues until full saturation occurs. 

In the case of an accretion disk, full saturation may not occur owing to the effects 
of turbulent viscosity in limiting the back-reaction. We therefore aim to find steady-state 
solutions in which the torque on the disk is nonzero. We also allow simultaneously for 
both aspects described above, by formulating a nonlinear problem in which the feedback is 
included self-consistently. In the process of the analysis, we will determine the size of the 
region over which the torque is injected into the disk. 

3. Analysis of the corotation region 

3.1. Basic equations 

We consider a three-dimensional gaseous disk having a barotropic equation of state, so 
that the pressure p depends only on the density p. We assume that the forcing potential and 
the properties of the unperturbed disk vary smoothly in radius across the corotation region. 
These assumptions may need to be reconsidered if the resonance is very close to the planet 
(i.e., at a distance comparable to the vertical scale- height H of the disk), or close to a sharp 
disk edge (where the density varies on a radial scale comparable to H) . However, the torque 
cut-off (GT79) limits the effectiveness of resonances that are closer than approximately H, 
and a gap in the disk about the planet's orbit may eliminate such close resonances in any 
case. We also consider the planet to execute a fixed orbit, and ignore any effects of migration. 
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The disk is governed by the equation of mass conservation, 

(d t + u ■ V)p = -pV • u, 

and the equation of motion, 

(dt + u- V)« = -V($ + /i) + -V • T. 

P 



(2) 



Here u is the velocity, $ the gravitational potential, /i = J p^ 1 dp the enthalpy, and T the 
stress tensor, which represents a viscous or turbulent stress. For a Navier-Stokes model we 
have 

T = /i[Vu+(VM) T ]+(/i b -f/i)(V-u)l, (3) 

where p is the shear viscosity and p^ the bulk viscosity. We assume that, as is conventional 
in accretion disk theory, the vertically integrated shear viscosity uT, may be expressed as 
a function of the radius r and the surface density S. The logarithmic derivative = 
d m(VE)/9 In E will feature in the analysis below. 

The unperturbed disk is a steady, axisymmetric solution {p^ u \ u^ u \ h^ u \ T^} of these 
equations in the presence of the steady, axisymmetric potential $( u ) of the central mass. Let 
the potential be perturbed by the addition of a nonaxisymmetric external potential $' that 
rotates rigidly with angular pattern speed fl p . The solution in the presence of the perturbed 
potential $ (p) = $ (u) + $' is {p (p) , u (p) , /i (p) , T (p) }, and may be assumed to be steady in a 
frame of reference that rotates at Jl p . 

We write the nonlinear equations for the Eulerian perturbations etc., 
in cylindrical polar coordinates (r, 0, z) in the rotating frame, using (u, v, w) to denote the 
cylindrical velocity components, and omitting the superscript (u): 



Dp' + (u'd r + w'd z )p = -p' 



-d r (ru) + d z w 



(P + p') 



-dJru') + -d^v' + d z vJ 
r r 



,/2 



Here 



Du + (u'd r + w'd z )u - 2Qv' - — = -d r (& + h!) + 

±D(rv') + 25m' = ~d+(& + h') + -->, 
Dw' + (u'd r + w'd z )w = -d g (& + h') + ---. 



D = (u + u')d r + ( Q - Q p + — ) + (w + w')d z , 



(4) 

(5) 

(6) 
(7) 

(8) 
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and the dots denote viscous terms that we do not write out owing to their complexity. 
Q = v/r is the unperturbed angular velocity in the inertial frame, which depends only on r 
in a barotropic disk, and 

2B = -^-(r 2 Q) (9) 

r dr 

is the vertical vorticity, or twice Oort's second parameter. The epi cyclic frequency k is given 
by k 2 = 40B. Also 

h'= V jp' + 0(p' 2 ), (10) 
where v 2 = dp/ dp is the square of the sound speed. 



3.2. Expansion about corotation 

Let r c be the corotation radius, defined by Q(r c ) = f2 p . There are three characteristic 
length-scales that are potentially relevant for determining the radial extent of the corotation 
region, 

/ \ 1/3 / lT , \ !/ 2 



5c k ' 5v {-mdQ/dr) ' {-K 2 dlnn/dkir) ' (U) 

Here c is an appropriate average of the sound speed v s , while m and ^ are the azimuthal 
wavenumber and amplitude of the forcing potential We assume that dVL/dr < 0. In a 
linear, inviscid theory (GT79) the scale is set by the evanescent wavelength 5 C of the two- 
dimensional mode (density wave) near corotation. However, the solution involves a cusp, 
which may be expected to be smoothed over a characteristic length 5 U in the presence of 
an effective viscosity. On the other hand, a ballistic treatment of the corotation region 
suggests that an annulus of librating trajectories is formed, of half-width 2^/25$, (Goldreich 
& Tremaine 1981). 

In order to allow all these effects to compete, we adopt scalings such that these three 
lengths are formally of the same order of magnitude. (Subsequently, we will take limits 
in which one or other of the quantities 5 can be distinguished as the relevant scale.) The 
characteristic width of the corotation region is then H, the semi-thickness of the disk, as 
was also found by Lubow (1990). This allows us to resolve the radial structure of the region 
in the same way as for the vertical structure, by introducing scaled coordinates 

e = ^, C=f, (12) 

where e 1 is a typical value of the angular semi-thickness H/r of the disk. In units such 
that r c = 1, the corotation region is found where £ and ( are of order unity. 
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The central potential may be expanded in the corotation region in a double Taylor 
series, 

$ = ($00 + e$ i£ + •••) + e 2 K 2 ( $ 2o + e$2i£ + •••) + •••■ (13) 

The perturbing potential may be expanded similarly, but is assumed formally to be smaller 
by a factor 0(e 2 ): 



$'a 2 = + e& 01 m +•••]+ ^ 2 K 2 [ $, 2o(0) + ^' 21 m + 



+ 



(14) 



This assumes that both potentials are symmetric about z = and vary smoothly in r in the 

neighborhood of r c . The scaling of $' will turn out to provide a critical level of nonlinearity 
in the solution. 

Properties of the unperturbed disk generally have non-trivial structure in ( but may be 
expanded in Taylor series in r, 

P = Po(C) + ePi(C)e + ---, 

u = e 3 MC) + •••], 

w = e 4 [w (C) + ■■■], 

Q = tt + efii£ H , 

B = B + eB^ + ---, 

v s = e[w s0 (C)H ], 

I 1 = e 3 [Po(0 + • • •] , 

Ph = e 3 [/i b0 (C) + • • •] • (15) 

Note that f2 = fi p by definition of the corotation radius. The scaling of \x implies that the 
dimensionless viscosity parameter, a = pQ/p, is formally 0(e), and will turn out to provide 
a critical level of dissipation in the solution. The scalings of the accretion flow, u and w, 
follow from this assumption. 

Perturbed variables generally have non-trivial structure in £ and (. We have found the 
required scalings to be 



P 

u 

V 

w 

h 
/i 

p'h 



= e[Po(£, 0,0 + •••], 

= e 2 K(0) + £«;(£, 0,o + •• 

= e > o (£,0) + e ^(£,0,c) + . 

= e 3 K(£,0,C) + •••], 

= e 3 [h' (C,<f>) + eh' 1 (C,<f>,0 + 

= e 4 K,(£,0,O + •••], 

= 6 4 [/i bo (£,0,O + ••■]• 



(16) 
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Equations (4), (5) and (6) at leading order, 0(e 2 ), then imply 

(«(,<% + Q^d^p'o + u' pi + w' d c p = -p ( — + d ( u[ + — <9^ + d C w 'oj > ( 17 ) 

-2Q ^ = -$' 01 - ^/i'o, (18) 

2£? < = (19) 
Also required is equation (6) at the next order, 0(e 3 ), 

(u' d ( + n 1 ^d <l ,)v , + 2B oU ' 1 + 2B 1 ^ = (^<t> i - -^)-Vo + -ft+— 

r c \ r c / r c PO PO 

(20) 

The last two terms are the only viscous terms to enter our analysis. Finally, equation (7) at 
leading order, 0(e 3 ), requires only a hydrostatic balance, 

= -$' 20 C - d c h' v (21) 

The enthalpy and density perturbations are related by equation (10) at leading order, 0(e 3 ), 

K = ^Vo- (22) 
Po 



3.3. Derivation of the governing equation 

We eliminate p' , u[, v' , and w' Q (by integrating eq. [17] with respect to ( over the full 
vertical extent of the disk) to obtain a single equation for h' , 

OQ / B \ 

(Ufa + fti£fy)(«jjJ 3 - hdl)h' + (hdl + 2r c Q Q 1 h)dlh' = —(9^) ( 7 2 - -±h , (23) 

involving the integrals 




(24) 

Here m' is given directly by equation (19). We recognize I\ = S and I 2 = Si as coefficients 
in the Taylor expansion of the unperturbed surface density S(r). Defining the mean sound 
speed c, the mean kinematic viscosity v, and the logarithmic viscosity derivative _D M in 
physical terms by 
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where 



we find I 3 = Eo/cq, 7 4 = z/ £ , and ^5 = -D^o£o/ c o- 

When the e-scalings and subscripts are omitted, the final equation in physical terms 
reads 

u' = (27) 

Here, all quantities except £, $'(0)5 M '(0)> an d 0) are to be understood as constants 
evaluated at r = r c . 

Equation (26) describes, in a reduced but systematically obtained manner, how the 
enthalpy perturbation is forced by the external potential in the corotation region. If the 
advective term involving v! were omitted, this would be a linear inhomogeneous equation, 
and the forcing would be proportional to the gradient of the vortensity H/B. The u' term 
derives from a product of perturbed quantities and is in this sense a nonlinear effect. An 
alternative way to write equation (26) is 

(28) 

where E' = E/i'/c 2 and IB' = dgv 1 are the leading-order perturbations of the surface density 
and the vertical vorticity. This shows that the nonlinear effect can be understood either as 
an additional advection of the perturbation or as the feedback of the perturbation on the 
vortensity gradient of the disk. 



3.4. Reduction to dimensionless form 

We assume that the external potential contains a single Fourier component, $' = 
^cosm^, and rewrite the governing equation (26) in a dimensionless form by means of 
the transformations 

f = 4 V( f ,« = in (§). (29) 

Thus 

(-a sin 9 8 x + xd e )(l - d 2 x )f - v(d 2 x - b)8 2 J = asintf, (30) 



9 



where 



a 



= 2 



V dlnQJ c 2 ' 



b = - 




mdVt/dr ) c 3 



(31) 



The problem therefore depends on only three dimensionless parameters: the forcing ampli- 
tude a, the viscosity derivative b, and the dimensionless viscosity v. For a Keplerian disk, 
b = 3-D^. In terms of the characteristic length-scales defined in equation (11), v = (S u /S c ) 3 . 

Equation (30) is to be solved on— oo<:r<oo,0<#<27r subject to a periodic 
boundary condition in 9, which implies that the physical solution has the same periodicity 
in as the forcing potential. The solution / should be bounded as \x\ — > oo in order to be 
a valid localized solution in the corotation region. 

We note immediately that one non-localized solution of equation (30) is / = — x + 
constant. This represents an axisymmetric perturbation that exactly cancels the vortensity 
gradient of the unperturbed disk, so nullifying the corotation resonance. Such a solution 
should be excluded, as it amounts to redefining the original problem. 



Equation (30) contains multiple x- derivatives but involves only one power of x explicitly. 
Like Airy's equation, it is amenable to Fourier analysis in x. Let 



be the Fourier transform of /, which exists provided that / decays sufficiently rapidly as 
\x\ — * oo. (In fact, by allowing for F to be a generalized function, we may permit / to grow 
at most algebraically in x.) The equation transforms to 



3.5. Solution by Fourier analysis 




(32) 




(33) 



where 




(34) 
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Equation (33) could be solved numerically as a partial differential equation, but is amenable 
to further Fourier analysis in 9. Let 

oo 

G(k,9)= G n (k)e me , (35) 

n=— oo 

then we obtain the system of ordinary differential equations 

+ ^ (irS) Gn ~ T (Gn+1 " Gn - l] = {<1 " Vl) W " (36) 

Note that the symmetry property 

G B (-fc) = - [G- n (k)]* (37) 

follows from the fact that the enthalpy perturbation h' is real. As a boundary condition, 
G n {k) is required to tend to zero as \k\ — > oo, in order that it be the Fourier transform of a 
continuous function. 

The total tidal torque exerted on the disk is 

/oo /»2"7r /»oo 
/ / (P + p')d^'rdrd(f)dz. (38) 
oo JO JO 

When expressed in terms of the scaled variables introduced in Section 3.4, the torque exerted 
in the corotation region is 



mrc 2 ^ d /£ 



2w roo 



where the prefactor is understood to be evaluated at r = r c . In terms of the Fourier 
representation, this further simplifies to 

T c = t c T GT , (40) 

where 

i c = Gi(0)-G_i(0) (41) 

is dimensionless, and 

mvr 2 ^ 2 d /E\ 

TGT= 2^V^)^UJ ( } 

is the torque formula given by GT79. Although equation (36) requires G±i(h) to have 
unit discontinuities at k — 0, the dimensionless torque t c is nevertheless well defined as the 
quantity G\{k) — G-\{k) is continuous at k — 0. 
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3.6. A simplifying assumption 

A considerable simplification in equation (36) is obtained when 6=1. This occurs in 
a Keplerian disk when conditions in the corotation region are such that z/X oc X 1 / 3 locally. 
In the absence of any detailed knowledge of the effective viscosity, we adopt this convenient 
assumption, anticipating that our results will not depend sensitively on it. 

It is then natural to rescale the wavenumber to 

k = u 1/3 k, (43) 

so that 

n d S^L + p Gn _ p ~ k ( Gn+1 _ G n -i) = (5 nA - <5 n ,_i) 5(k), (44) 
dk 

where 

p = \aV~^ (45) 
is the only remaining parameter of the problem. In physical terms, 

V ~ \-K 2 d\nn/d\nr) \ v ) ( ' 

measures the strength of the potential relative to the viscosity (i.e., the nonlinearity relative 
to the dissipation). In terms of the three characteristic length-scales defined in equation 

(ii), P = (<W^) 2 . 

We may assume without loss of generality that p > 0. Our aim now is to evaluate the 
function t c (p) that determines how the formula of GT79 is modified. 



3.7. Linear theory 

In linear theory the coupling parameter p is set to zero and only modes n — ±1 are 
excited. The solution is 

G ±1 (k) = ±H(±k) exp( T §P), (47) 

where H is the Heaviside function. The dimensionless torque (eq. [41]) is then t c — l. This 
shows that the torque formula of GT79, which was derived for a two-dimensional, inviscid 
disk, also applies to a three-dimensional, barotropic and viscous disk, provided that the 
linearization is valid (i.e., p <^ 1). 
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3.8. Feedback 

The axisymmetric component Go of the solution has a special role. The case n = in 
equation (44) is not a differential equation but states that 

~k 2 G =p~k{G 1 -G- 1 ). (48) 

Now G also appears in the equations for G±\ in the form pkG . We solve equation (48) 
within the space of generalized functions to obtain 

kG =p(G 1 -G- 1 ) + C 1 6(k), (49) 

where C\ is an arbitrary constant. Equation (44) with n = ±1 then implies 



dk 
dG_i 



+ k I G l - P kG 2 + p'id - G_i) = (1 - C lP )5(k), (50) 
- PGU - pkG- 2 + p 2 (G 1 - GLi) = (1 - Cip)<J(fc). (51) 



eta 

The solution of equation (49) is of the form 

Go = G r cg (k) + & - CKJ'(fc) + ^G 2 5(fc), (52) 

where GQ Cg (A;) is regular at k — 0, and G 2 is a second arbitrary constant. Consider the inverse 
Fourier transforms 

/OO /"OO /"OO 

A;- 1 e ifc:c dA; = i7rsgn(a;), -/ 5\k) e ikx dk = ix, / i 5(k) e lkx dk = i, (53) 
-OO J —OO J —OO 

the first of which is a principal value. By allowing Ci 7^ we would admit a perturbation 
that redefined the basic vortensity gradient of the disk across the whole corotation region. 
This should be excluded, as discussed in Section 3.4. We may also set C2 = without loss 
of generality as it simply adds a constant to the basic surface density. 

The solution generated by the k~ l singularity in Go corresponds to a gradual, axisym- 
metric step in the enthalpy perturbation across the corotation region, in the sense that f(x, 6) 
tends to well defined limits as x — > ±00, but these limiting values are unequal. This occurs 
because the tidal torque causes a redistribution of the surface density of the disk. Consider 
the standard evolutionary equations for an accretion disk in which the corotation torque is 
added as an infinitely localized source, 

f + ;|(rav) = 0, (54) 
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d ( 2 . 19/ gfif^ T c 



E^-(r 2 ft) = -— z/Er" 3 — + — ^- 5(r - r c ), (55) 
ar r or \ ar J 2irr c 

where v r is the mean radial velocity. In a steady state these integrate to 

— 2nrT,v r — M — constant, (56) 

Mr 2 tt + 27rz/Er 3 ^ + T c H(r - r c ) = constant. (57) 
ar 

The step in z/E across the resonant region associated with the k~ x singularity in G can be 
shown exactly to balance the corotation torque T c in this equation. Note that this simplified 
description does not treat the internal structure of the corotation region, but is intended to 
show how the torque affects the disk on a larger scale. 



3.9. Approximate solutions for small and large p 

Approximate solutions can be obtained by neglecting G2 and G_ 2 in equations (50) and 
(51). This amounts to a severe Fourier truncation of the problem, in which we allow for 
the solution to feed back on the axisymmetric part of the disk but neglect the excitation of 
higher harmonics. Nevertheless, this method is found to give results in close agreement with 
the more accurate numerical solutions described in the next section. 

Defining G± — G\ ± G_i and taking the sum and difference of equations (50) and (51), 
we obtain 

AC 

" v ' 4 " + (P + 2p 2 )GL = 26(h), (58) 



dk 



and so 



dk 

d ( 1 dG 



■ k-G. 0. (59) 



+ (r + 2p 2 )G- = 28(h). (60) 



dk \k 2 dk 

For p 2 1 the solution may be expanded in powers of p 2 , 

G_ = G {0) +p 2 G m +0(p 4 ). (61) 

Using the fact that the bounded solution of the equation 

d 2 y 



is 

y(x) 



, r2 +y = 2f(x) (62) 

/oo 
f(x') e-l*-*'l dx' : (63) 
-00 
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we find, at successive orders in p 2 , 

G { °\k) =exp(-i|jfe 3 |), (64) 

POO 

G { l\k) = - / exp(-i|fc ,3 | - i|jfe 3 - k f3 \) dk'. (65) 



(66) 



The dimensionless torque is then 

2/3 



t c = G_(0) = !-(|) r Q) + °(p 4 ) ~ 1 " 2 - 044 ^ 2 - 



This asymptotic form can be shown to be correct even when the coupling to higher har- 
monics is taken into account. In this limit of small p, only wavenumbers \k\ < 1 contribute 
significantly to the solution. 

For p 2 > 1 we may neglect k relative to 2p 2 in equation (60). The solution is then 

G-ik) « ^p~^\kf 2 K m (2-^pk 2 ) , (67) 

where K is the modified Bessel function of the second kind, and the dimensionless torque is 

r (1\ 

t c = G_(0) w — ^-2 1 /V 3/2 ~ 0.4019 p~ 3/2 . (68) 

Since p oc z/ -2 / 3 , this implies that the torque is proportional to the viscosity when the 
viscosity is small. This is to be expected on general grounds, but to neglect the coupling 
to higher harmonics is difficult to justify formally when p is large. The above form of G_ 
indicates that only wavenumbers \k\ < p^ 1 ^ 2 contribute significantly. 



3.10. Numerical solution 

For a numerical solution, we truncate the system at some order iV and set G n to zero for 
\n\ > N. We then solve equations (50) and (51) together with equation (44) for 2 < \n\ < N. 
Now the symmetry property (37) allows us to consider k > only, and the jump conditions 
at k — imply 

G n (0+) + G_ n (0+) = 6 n , u n>0. (69) 

We integrate the equations from k = 0+ to a finite value k = k max . A 'particular solution' 
Gn^ satisfying the inhomogeneous boundary condition (69) is generated by starting from 
the initial condition G 1 — 1 (G n = otherwise), and iV 'complementary functions' Gn'^ 
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by starting from initial conditions G q = —G- q = 1 (G n = otherwise) for each q = 
1, . . . , N. The amplitudes of the complementary functions appearing in the desired solution 
are determined by requiring that G„(A; max ) = for — N < n < — 1. This simulates the 
requirement that G n should tend to zero as \k\ — > oo. For n < 0, G n has a tendency to 
grow as exp(|A; 3 /|n|) as k — > oo, and our boundary condition is designed to eliminate such 
growing solutions. For n > 0, G n tends to decay naturally as k — > oo. 

In Figure 1 we plot the dimensionless torque t c determined from the numerical solution 
as a function of the coupling parameter p. The solution converges rapidly as iV and k max 
are increased, and good agreement is found with the approximate solutions (66) and (68) for 
small and large p. 



3.11. Interpretation 



We now attempt to give a real-space interpretation of these results. Multiplying our 
governing equation (26) by — £r/(2f2) and integrating over 0, we obtain 



a? 



vr 
2fi 







^ dVL D,. . , 
2rn- £ ) ti 



d(f> 



(d^')d ( h' #, (70) 



where the very first term has been taken over to the right-hand side. The first integral with 
respect to £ is 



plus a possible constant of integration. 



} ^="§/ (71) 



Equation (71) expresses the conservation of angular momentum in a steady state, al- 
though the terms that balance in the unperturbed disk do not appear. The left-hand side 
is the divergence of the angular momentum flux, and the right-hand side is the tidal torque 
per unit radius. To see this, the first term in the flux can be identified as the integrated 
Reynolds stress 

r J J pu'v' d(pdz, (72) 

using equations (18) and (19). The remaining terms in the flux are viscous stresses. The 
right-hand side of equation (71) can also be written as 



p'dff,®' d(f) dz, 



(73) 



which is the tidal torque per unit radius. 
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In the Fourier representation, when suitably normalized, equation (71) reads 

^(F 1 - F_x) - vk{k 2 + b)F = —(F, - F_ x ), (74) 

plus a possible multiple of <5(/c). The same equation could be written for G instead of F, 
and is then equivalent to our equation (36) in the case n — 0. (The possible additional term 
proportional to S(k) is equivalent to the C\ term that arises in passing from equation (48) to 
(49), and can be excluded on the same grounds.) This shows how equation (36) partitions 
into terms associated with Reynolds stress, viscous torque and tidal torque. 

In Section 3.9 we worked with the quantity G_ oc (1 + k 2 )(Fi — F_i). In these terms, 
the Fourier transforms of the Reynolds, viscous and tidal terms are proportional to 

(ttp) g - w - G - {k) ' (ttp) g - w - (75) 

respectively. When p <C 1 we found that G_(k) is cut off for wavenumbers \k\ > 1, i.e. 
physical wavenumbers \k/8 c \ > 1/S V . It follows that the characteristic scale of the viscous 
torque is the viscous scale 5 U defined in equation (11). If the viscosity is small, so that 
5 U <C S c , then the scale of the Reynolds stress is also 5 U . However, the scale of the tidal 
torque is 5 C , because its Fourier transform is cut off for physical wavenumbers \k/5 c \ > l/S c 
by the factor (1 + k 2 )' 1 . 

When p > 1 we found that G-(k) is cut off for wavenumbers \k\ > p _1//2 , i.e. physical 
wavenumbers \k/S c \ > l/5y. If the forcing amplitude is small, so that 5^ <C S c , then the 
scale of the Reynolds stress is also 5y, but the scale of the tidal torque is again S c . 

To summarize, if 8 C is the largest of the three scales, then the tidal torque is spread over 
a region of scale 5 C , while the Reynolds and viscous stresses have a smaller scale max(<5„, 5^). 
This is consistent with the impression obtained from the original theory (GT79), that the 
tidal torque creates a disturbance in the disk on the scale 5 C , which subsequently transfers its 
angular momentum to the background disk, through dissipation or nonlinearity, on a finer 
scale. 



4. Application to eccentric resonances 

The reduction of the torque exchanged between the perturber and the disk has con- 
sequences for the eccentricity evolution of young planets orbiting in a protoplanetary disk. 
We consider the first-order eccentric corotation resonances associated with a planet of mass 
ratio q = M p /M* executing a orbit of eccentricity e within a Keplerian disk. 



-17- 



GT80 provide expressions for the locations of these resonances and the forcing potentials. 
Inner and outer resonances occur at radii 

/ \ 2/3 

where a p is the semi-major axis of the planet's orbit. The amplitude of the forcing potential 
may be written in the form 

0.8020 C±meq, (77) 



r 2 tt 2 

where is a correction factor that tends to unity for large m. Accordingly, the saturation 
parameter of our analysis is 

/rr\ -4/3 

p = 0.7006 C±m 5/3 eqa- 2/3 - , (78) 

where we write v = acH and H = c/Q. When p — 1, i.e., when 

1.427 , 2/ , //A 4/3 . . 

e = eo.es = -^m-^q- 1 ^ 3 {-J , (79) 

63% saturation occurs because t c 0.37. For 5% saturation, only p m 0.16 is required. 

Table 1 contains results for the lowest-order outer eccentric resonances (l — m — 1 
in the notation of GT80) for a disk with H/r = 0.05 and a planet of mass ratio q = 
0.001. The first column labels the azimuthal wavenumber of the potential component being 
considered. The second column gives the location of the resonance in units of a p . The third 
column contains the correction factor defined in equation (77). The fourth column gives a 
measure of the relative importance of the resonances: neglecting the torque cut-off and any 
differences in the disk parameters between different resonant radii, the unsaturated torque 
scales oc m\l/ 2 oc C^m 3 . The fifth column, labeled eo.63, is the eccentricity required to reduce 
the eccentric corotation torque 63% below its unsaturated value for a turbulent viscosity 
parameter a = 0.004. The sixth column, eo.os, is similar to the fifth column, but gives the 
eccentricity needed to reduce the torque 5% below its unsaturated level. The final column, 
a .63, is the value of a corresponding to 63% saturation when the eccentricity is e = 0.1. 
Table 2 gives the same data for the inner resonances (I — m + 1). 

Resonances of high order occur very close to the planet's orbit, and are present only 
if the planet does not clear a substantial gap in the disk. An analysis of such resonances 
requires a consideration of a number of effects that we have neglected. When m becomes 
comparable to r/H, additional terms in the dynamical equations must be included, and the 
torque is reduced from the standard value (GT79). The proximity of the resonance to the 
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planet also means that the perturbing potential varies significantly with z within the disk, 
and the relevant quantity is the density-weighted vertical average of the potential, rather 
than the mid-plane value (Ward 1986, 1989). These two effects apply equally to Lindblad 
and corotation resonances, and so do not alter the 5% balance between torques that excite 
and damp eccentricity. A further effect is that neighboring resonances of sufficiently high 
order overlap one another, which may alter the resonant torques. 

5. Summary and discussion 

We have determined the torque exerted in a steady state by an external potential on 
a three-dimensional gaseous disk at a non-coorbital corotation resonance. Our model ac- 
counts for the feedback of the torque on the surface density and vorticity in the corotation 
region, and assumes that the disk has a barotropic equation of state and a nonzero effective 
viscosity. The torque formula of Goldreich & Tremaine (1979) (eq. [42]) must be modified 
by a reduction factor t c (p), plotted in Figure 1, which quantifies the extent to which the 
resonance is saturated. This factor depends essentially on a single dimensionless parameter 
p, defined in equation (46), which measures the strength of the potential \I/ relative to the 
viscosity v (i.e., the nonlinearity relative to the dissipation). 

In Section 3.10, we determined that the characteristic radial width of the resonance in 
the limit of large p (low viscosity) is 5 ~ V^/Q. (This refers to the scale of Reynolds and 
viscous stresses; the tidal torque is spread over a region of characteristic width c/k.) We 
found that the torque is reduced by a factor ~ p~ 3 / 2 , which can be regarded as the ratio of 
the viscous diffusion rate ~ u/5 2 over the libration region to the libration rate ~ m\dQ/dr\5. 
In this regime, our results are broadly consistent with the saturation model of Ward (1992), 
which involves the same ratio of rates, but applied to the coorbital region. Earlier, Goldreich 
& Tremaine (1981) had argued that, in a disk of collisional particles, saturation would occur 
when (in our notation) p ^> 1. 

Goldreich & Sari (2002) have recently developed an evolutionary model for planetary 
eccentricity in which eccentricity growth occurs through a finite-amplitude instability. For 
infinitesimal eccentricity, the damping caused by corotation resonances just overcomes the 
eccentricity growth due to Lindblad resonances. Above a critical level of eccentricity, how- 
ever, the corotation resonances become sufficiently saturated that growth occurs. Using our 
evaluation of the saturation function t c (p), Goldreich & Sari (2002) were able to determine 
the critical value of eccentricity required within the context of a specific disk model. The 
results in Section 4 and Table 1 show that for a Jupiter-mass planet orbiting within a typical 
protoplanetary disk, eccentricities of a few percent are adequate to saturate all first-order 
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eccentric corotation resonances except those of the lowest m- values (m < 4). To achieve a 
5% reduction in torque at such resonances, as may be sufficient to change the balance from 
eccentricity damping to growth (GT80), requires only eccentricities of 1% or less. 

Exactly which resonances contribute most to eccentricity evolution depends on the ex- 
tent of the gap cleared by the planet, which in turn depends on the mass ratio and the 
properties of the disk. As seen in Table 1, the largest corotation torques are those associated 
with the highest m-values, and they are the easiest to saturate. A lower-mass planet that 
opens a smaller gap may excite many resonances so that the dominant torque comes from 
approximately the cut-off m-value, of order r/H (GT80). In that case, 63% saturation is 
achieved optimistically at the torque cut-off when e > 1.4 g _1 o; 2//3 (iJ/r) 3 . Note that, al- 
though the large value of m is advantageous for saturation, the low mass of the planet is 
unfavorable (see eq. [79]). On the other hand, if the planet mass is small enough so that 
there is no gap in the disk, then the coorbital resonances cause eccentricity decay (Ward 
1988; Artymowicz 1993). 

We have neglected a number of potential complications. Numerical simulations of a 
Jupiter-mass planet in a circular orbit with the disk parameters adopted in Section 4 suggest 
that the disk edge is not sharp (i.e., its radial extent is much greater than H; see Figure 1 
of Lubow, Seibert, & Artymowicz 1999). It is evident, however, that within the radial 
extent of the planet's Roche lobe, material is captured by the planet and the resonant effects 
considered here do not play a role. It is not clear which eccentric corotation resonances are 
excited in the presence of an eccentric planet. The complications are due at least in part to 
nonlinear effects other than those considered here. For example, shocks associated with the 
planet's wake cause a non-closure of streamlines in the vicinity of the disk edge, leading to 
a drift of material through the corotation regions. In addition, material within the Roche 
lobe of the planet exerts torques that may contribute to the eccentricity balance. 

As usual in accretion disk theory, our knowledge of the effective viscosity of the disk 
is limited. In Section 3.6 we adopted a convenient assumption in order to make analytical 
progress. Although we have succeeded in analyzing a three-dimensional disk, thereby gener- 
alizing existing theories, we assumed that the disk was barotropic. The effects of buoyancy 
or baroclinicity on the corotation region remain uninvestigated. 

In recent work we were able to verify analytical theories of the torques exerted at 
Lindblad and vertical resonances through direct numerical simulations (Bate et al. 2002). 
It would be valuable to conduct simulations of a non-coorbital corotation resonance to test 
the findings of the present paper. 

After eccentric corotation resonances saturate at small eccentricity, the eccentricity 
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growth may be limited at intermediate values. The limitation may be due to the overlap 
of resonances or alternatively through the excitation of higher order Lindblad resonances, 
some of which cause eccentricity damping. A contribution to eccentricity damping occurs 
at an eccentric inner (outer) Lindblad resonance that lies outside (inside) the planet's orbit. 
SPH simulations of eccentric-orbit binary stars suggest that little eccentricity growth via 
resonances occurs for eccentricities in the range of 0.5 — 0.7 or higher (Lubow & Artymowicz 
1993), and similar limits may occur for planets. 

In conclusion, we have found a simple quantitative measure of the saturation of a coro- 
tation resonance in a gaseous disk. This analysis suggests that planets may plausibly ex- 
perience a net growth of eccentricity through their interaction with the disk in a variety of 
circumstances, provided that the eccentricity is not extremely small to begin with. 

We are grateful to Peter Goldreich and Re'em Sari for allowing us to coordinate our 
presentation of results with theirs prior to publication, and also for correcting some inac- 
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initiated at the Aspen Center for Physics, and we are grateful for their hospitality. GIO 
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Table 1. Saturation parameters for outer corotation resonances 



m 


r/a p 


c~ 


C> 3 


e 0.63 


eo.05 


«0.63 


2 


1.5874 


0.7422 


4.407 


0.2812 


0.0453 


0.00085 


3 


1.3104 


0.8418 


19.13 


0.1261 


0.0203 


0.00282 


4 


1.2114 


0.8854 


50.18 


0.0742 


0.0120 


0.00625 


5 


1.1604 


0.9101 


103.5 


0.0498 


0.0080 


0.01139 


6 


1.1292 


0.9261 


185.2 


0.0361 


0.0058 


0.01843 


7 


1.1082 


0.9372 


301.3 


0.0276 


0.0044 


0.02759 


8 


1.0931 


0.9454 


457.6 


0.0219 


0.0035 


0.03903 


9 


1.0817 


0.9517 


660.3 


0.0179 


0.0029 


0.05292 


10 


1.0728 


0.9567 


915.3 


0.0149 


0.0024 


0.06941 



Table 2. Saturation parameters for inner corotation resonances 



m 


r/a p 


Cm 


C> 3 


eo.63 


eo.05 


«0.63 


1 


0.6300 


0.3365 


0.113 


1.9688 


0.3170 


0.00005 


2 


0.7631 


1.1819 


11.17 


0.1766 


0.0284 


0.00170 


3 


0.8255 


1.1265 


34.26 


0.0942 


0.0152 


0.00437 


4 


0.8618 


1.0970 


77.02 


0.0599 


0.0096 


0.00862 


5 


0.8855 


1.0787 


145.5 


0.0420 


0.0068 


0.01469 


6 


0.9023 


1.0663 


245.6 


0.0314 


0.0050 


0.02277 


7 


0.9148 


1.0572 


383.4 


0.0245 


0.0039 


0.03305 


8 


0.9245 


1.0503 


564.8 


0.0197 


0.0032 


0.04570 


9 


0.9322 


1.0449 


795.9 


0.0163 


0.0026 


0.06088 


10 


0.9384 


1.0406 


1083. 


0.0137 


0.0022 


0.07873 
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Fig. 1. — The dimensionless torque versus the coupling parameter. The curves calculated 
for truncation orders N = 1, ... ,5 and /c max — 3 are plotted as solid lines; the curve for 
N — 1 lies slightly below the others, while those for N > 1 are indistinguishable by eye. 
The dotted line shows the small-p approximation (66) and the dashed line shows the large-p 
approximation (68), to which the curve for N = 1 asymptotes. 



